################################################################################
#this is a script to create weather trace data for specific locations

#script was written by Jeremy VanDerWal (jjvanderwal@gmail.com), 6 November 2008
#and modified by Luke Shoo 23rd June 2009
#and modified by Collin Storlie 1st November 2010

################################################################################
#load necessary libraries

#list the libraries needed
necessary=c("adehabitat","SDMTools")

#check if library is installed
installed = necessary %in% installed.packages()

#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)

#load the libraries
for (lib in necessary) library(lib,character.only=T)

################################################################################

#give the necessary directories
tmax.dir = "/home1/99/jc152199/In/AWAP/"

#define the output folder
outfolder = "/home1/99/jc152199/Out/AWAP/"

#define the starting date
start.day = as.Date("2006-11-11","%Y-%m-%d")

#define the last date
last.day = as.Date("2009-06-01","%Y-%m-%d")

#define your sites ... ensure locations are not on the max and min of longitude or latitude
locationdata=read.csv("/home1/99/jc152199/In/Logger Locs for Downscale.csv")
X=locationdata$longdec
Y=locationdata$latdec
pnts = cbind(X,Y)

#prepare outdata
#fill the julian day column with all possible days from the start to finish date
outdata = data.frame(Julian=0:julian(last.day,origin=start.day),date=NA, AU10A2=NA, AU1A3=NA, AU2A2=NA, AU4A2=NA, AU6A2=NA, AU8A2=NA, BK10A2=NA, BK12A2=NA, BK14A2=NA, BK16A2=NA, CU10A2=NA, CU12A2=NA, CU1A1=NA, CU2A2=NA, CU4A2=NA, CU6A2=NA, CU8A2=NA, SU10A2=NA, SU10B2=NA, SU350A2=NA, SU6A2=NA, SU8A2=NA, WU9A2=NA, WU11A2=NA, WU13A2=NA, S31011=NA, S31034=NA, S31052=NA, S31108=NA, S31193=NA, S31209=NA, S31210=NA, S32004=NA, S32025=NA, S32037=NA, S32040=NA, S32078=NA, S32141=NA)

#convert the julian day column to a date column
outdata$date = format(as.Date(outdata$Julian,origin=start.day),"%Y%m%d")

################################################################################

#get a list of years
years=2006:2009

#Cycle through the years, months and days
for (year in years) {
  
  #doy = 0 #track the number of days in the year
  for (month in 1:12) {
    
    #dom = 0 #track the number of days in the month
    for (day in 1:31) {
      
      #calculate the date
      tdate = format(as.Date(paste(year,month,day,sep="-"),"%Y-%m-%d"),"%Y%m%d")
      
      #if the date exists as a file... start the summary
      tfile = paste(tmax.dir,year,"/tmax.",tdate,tdate,".grid.gz",sep="")
      if (file.exists(tfile)){  cat(tfile,'\n')
      
        #unzip the file and import the asc file
        tasc = read.asc.gz(tfile)
        
          #ensure the layer has data in it
          if (length(which(!is.na(tasc)))>0) {

            #extract the location specific info
            info = extract.data(pnts,tasc)
    
              #COLLIN MODIFY - IMPORTANT, UPDATE THE NUMBER OF COLUMNS OF DATA YOU WANT ADDED [ie. 3:? (2 + no of sites)]  
              outdata[which(outdata$date==tdate),3:40] = info
          }
        }
      }
    }
  }

#need to write out the summary data

write.csv(x=outdata,file="/home1/99/jc152199/Out/AWAP/AWAP_out.csv", sep="", row.names=F)
write.csv(x=sample.out,file=paste(out.dir, "sample.intersect.csv", sep=""), row.names=F)

sub1=outdata[,1:200]
subdata=outdata[,201:399]
sublab=outdata[,1:2]
sub2=cbind(sublab,subdata) 

write.csv(x=sub1,file="C:/LPS_working/Wellbergen/Tmax_summary_Justin_sites_1950_to_presentA.csv")
write.csv(x=sub2,file="C:/LPS_working/Wellbergen/Tmax_summary_Justin_sites_1950_to_presentB.csv")